home *** CD-ROM | disk | FTP | other *** search
/ USGS: Oil & Gas Potential…National Wildlife Refuge / USGS - Oil & Gas Potential of the Arctic National Wildlife Refuge - Disc 2.iso / mac / MEcode / MEPU.for < prev    next >
Text File  |  1999-02-11  |  5KB  |  407 lines
  1. c   MEPU.for generates uncertainty estimates
  2.  
  3. c    at 95, 50, 5th fractiles and the corresponding size
  4.  
  5. c    distribution.  Note, distribution is average of size at
  6.  
  7. c    fractile +/- 10 observations on either size.
  8.  
  9. c    These uncertainty estimates are for recoverable oil
  10.  
  11. c    or NA gas for ANWR individual plays. Program could be easily
  12.  
  13. c    modified to get in-place estimates of uncertainty.
  14.  
  15. c
  16.  
  17. c    Written by Schuenemeyer 3/22/98
  18.  
  19. c
  20.  
  21. c   Input files:
  22.  
  23. c    Unit  8 - play file
  24.  
  25. c    Unit  9 - prospect file
  26.  
  27. c   Output files:
  28.  
  29. c    Unit 10 - uncertainty estimates
  30.  
  31. c     Note: output will be pasted into ANWR Summary & Distns
  32.  
  33. c      Excel spreadsheets
  34.  
  35. c
  36.  
  37. c   Subroutines required:
  38.  
  39. c     SizeClass, Buble, Bublei(included in this program)
  40.  
  41. c
  42.  
  43.       character*25 fnpl,fnpr,fndn
  44.  
  45.     character title*80
  46.  
  47.       dimension x(10000),id(10000),tmp(7),nar(30),naid(30),so(2,13)
  48.  
  49.     dimension nt(4),sm(2),perc(3)
  50.  
  51.     integer fsc
  52.  
  53.     data perc/0.05,.5,.95/
  54.  
  55.     call getdat(iyrx,imonx,idayx)
  56.  
  57.     marry=13
  58.  
  59.     nsper=10
  60.  
  61.       write(*,*)' Enter play file name'
  62.  
  63.     read(*,'(a25)')fnpl
  64.  
  65.       open(8,file=fnpl)
  66.  
  67.     write(*,*)' Enter prospect file name'
  68.  
  69.     read(*,'(a25)')fnpr
  70.  
  71.       open(9,file=fnpr)
  72.  
  73.     write(*,*)' Enter output uncertainty file name'
  74.  
  75.     read(*,'(a25)')fndn
  76.  
  77.       open(10,file=fndn)
  78.  
  79.     7    write(*,8)
  80.  
  81.     8 format(/' Enter total number of PLAYS & Oil/Gas(1or2)')
  82.  
  83.     read(*,*) npl,nog
  84.  
  85.     if(nog.lt.1.or.nog.gt.2) goto 7
  86.  
  87.     if(nog.eq.1) then 
  88.  
  89.      scdiv=1.
  90.  
  91.      idv=1
  92.  
  93.       else
  94.  
  95.      scdiv=6.
  96.  
  97.      idv=3
  98.  
  99.     end if
  100.  
  101.     knt=0
  102.  
  103. c  read oil or gas values from play file
  104.  
  105.     read(8,*)
  106.  
  107.    10    read(8,*,end=20)it,(tmp(i),i=1,7)
  108.  
  109.     knt=knt+1
  110.  
  111.     id(knt)=it
  112.  
  113.     x(knt)=tmp(idv)
  114.  
  115.     goto 10
  116.  
  117.    20    close(8)
  118.  
  119.       iexc=npl-knt
  120.  
  121.     if(iexc.lt.0) pause 1
  122.  
  123.     nbs=0
  124.  
  125. c  big percentile (fractile) loop
  126.  
  127.       do ki=1,3
  128.  
  129.       ndep=0
  130.  
  131.       do i=1,2
  132.  
  133.      sm(i)=0
  134.  
  135.      do j=1,marry
  136.  
  137.       so(i,j)=0.
  138.  
  139.      end do
  140.  
  141.     end do
  142.  
  143.     maxf=0
  144.  
  145.       if(nbs.eq.0) then
  146.  
  147. c   sort oil/gas play totals
  148.  
  149.      call buble(x,id,knt)
  150.  
  151.     nbs=1
  152.  
  153.     end if
  154.  
  155.       npid=perc(ki)*npl+.5
  156.  
  157. c   percentile is zero
  158.  
  159.     if(npid.le.iexc) goto 56
  160.  
  161. c   relate npid to actual data and interval bounds
  162.  
  163.     ndat=npid-iexc
  164.  
  165.     une=x(ndat)
  166.  
  167.      write(*,*) ' Uncert est',une
  168.  
  169. c   get size distribution
  170.  
  171.       nll=ndat-nsper
  172.  
  173.     if(nll.lt.1)nll=1
  174.  
  175.     nul=ndat+nsper
  176.  
  177.     if(nul.gt.knt)nll=knt
  178.  
  179.     ndif=nul-nll+1
  180.  
  181.     xndif=ndif
  182.  
  183. c   store id's
  184.  
  185.       do i=nll,nul
  186.  
  187.      ia=i-nll+1
  188.  
  189.      nar(ia)=id(i)
  190.  
  191.      naid(ia)=ia
  192.  
  193. c    write(*,*)ki,i,ia,nar(ia)
  194.  
  195.     end do
  196.  
  197.       call bublei(nar,naid,ndif)
  198.  
  199. c   nar contains the sorted number for the prospect file.
  200.  
  201. c   Get play data and summarize.
  202.  
  203. c
  204.  
  205. c   read headers from prospect file
  206.  
  207.     rewind 9
  208.  
  209.       read(9,'(a70)')title
  210.  
  211.     read(9,*)nvar,nobs
  212.  
  213.     do i=1,nvar
  214.  
  215.      read(9,*)
  216.  
  217.     end do
  218.  
  219.     nsr=0
  220.  
  221.       do i = 1,ndif
  222.  
  223.      ia=nar(i)
  224.  
  225.     if(nsr.eq.1) goto 30
  226.  
  227.    28    read(9,*,end=59) (nt(k),k=1,4),oilg
  228.  
  229.    30 nsr=0
  230.  
  231.     if(nt(1).lt.ia) goto 28
  232.  
  233.     if(nt(1).gt.ia) then
  234.  
  235.      nsr=1
  236.  
  237.      goto 33
  238.  
  239.     end if
  240.  
  241.     if(nt(3).eq.1) nt2=nt(2)
  242.  
  243.     if(nt(4).eq. nog) then 
  244.  
  245.     pet=oilg/scdiv
  246.  
  247.     call SizeClass(pet,maxf,marry,fsc)
  248.  
  249. c  accumulate
  250.  
  251.     ndep=ndep+1
  252.  
  253.       so(1,fsc)=so(1,fsc)+1
  254.  
  255.     so(2,fsc)=so(2,fsc)+oilg
  256.  
  257.     end if
  258.  
  259.     nt2=nt2-1
  260.  
  261.     if(nt2.gt.0)goto 28
  262.  
  263.    33    end do
  264.  
  265. c  summarize and write out results
  266.  
  267.     do j=1,marry
  268.  
  269.     do i=1,2
  270.  
  271.      so(i,j)=so(i,j)/xndif
  272.  
  273.      sm(i)=sm(i)+so(i,j)
  274.  
  275.     end do
  276.  
  277.     end do
  278.  
  279.    56    fra=(1.-perc(ki))*100.
  280.  
  281.       write(10,60)fra
  282.  
  283.    60 format(f15.0)
  284.  
  285.    61    do i=1,2
  286.  
  287.     write(10,65)sm(i),(so(i,j),j=1,marry)
  288.  
  289.    65 format(3x,14f12.3)
  290.  
  291.       end do
  292.  
  293.       write(*,*) ' Max size class',maxf
  294.  
  295.     end do
  296.  
  297.       close(10)
  298.  
  299.     stop
  300.  
  301.    59 pause 3
  302.  
  303.       END
  304.  
  305.  
  306.  
  307.     Subroutine SizeClass(size, maxf, marray, fsc)
  308.  
  309. c Computes size class in log based 2
  310.  
  311.       integer fsc
  312.  
  313.        fsc = Int(Log(size) / Log(2.)) - 1
  314.  
  315.        If (fsc .lt. 1) fsc = 1
  316.  
  317.        If (fsc .gt. maxf)  maxf = fsc
  318.  
  319.        If (fsc .gt. marray)  fsc = marray
  320.  
  321.     return
  322.  
  323.       End 
  324.  
  325.  
  326.  
  327.     SUBROUTINE BUBLE(X,ID,N)
  328.  
  329.       DIMENSION X(1),ID(1)
  330.  
  331.       KS=N
  332.  
  333.    15 KW=0
  334.  
  335.       DO 30 I=2,KS
  336.  
  337.       IF(X(I).GE.X(I-1)) GOTO 30
  338.  
  339.       TMP=X(I)
  340.  
  341.       X(I)=X(I-1)
  342.  
  343.       X(I-1)=TMP
  344.  
  345.       NTI=ID(I)
  346.  
  347.       ID(I)=ID(I-1)
  348.  
  349.       ID(I-1)=NTI
  350.  
  351.       KW=1
  352.  
  353.    30 CONTINUE
  354.  
  355.       IF(KW.EQ.0) RETURN
  356.  
  357.       KS=KS-1
  358.  
  359.       IF(KS.EQ.1) RETURN
  360.  
  361.       GOTO 15
  362.  
  363.       END
  364.  
  365.  
  366.  
  367.     SUBROUTINE BUBLEI(X,ID,N)
  368.  
  369.       integer x,tmp
  370.  
  371.       DIMENSION X(1),ID(1)
  372.  
  373.       KS=N
  374.  
  375.    15 KW=0
  376.  
  377.       DO 30 I=2,KS
  378.  
  379.       IF(X(I).GE.X(I-1)) GOTO 30
  380.  
  381.       TMP=X(I)
  382.  
  383.       X(I)=X(I-1)
  384.  
  385.       X(I-1)=TMP
  386.  
  387.       NTI=ID(I)
  388.  
  389.       ID(I)=ID(I-1)
  390.  
  391.       ID(I-1)=NTI
  392.  
  393.       KW=1
  394.  
  395.    30 CONTINUE
  396.  
  397.       IF(KW.EQ.0) RETURN
  398.  
  399.       KS=KS-1
  400.  
  401.       IF(KS.EQ.1) RETURN
  402.  
  403.       GOTO 15
  404.  
  405.       END
  406.  
  407.